Properties of Foreshocks and Aftershocks of the Non-Conservative SOC 
Olami-Feder-Christensen Model: Triggered or Critical Earthquakes? 
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Following Hergarten and Neugebauer ,l| who discovered aftershock and foreshock sequences in 
the Olami-Feder-Christensen (OFC) discrete block-spring earthquake model, we investigate to what 
degree the simple toppling mechanism of this model is sufficient to account for the properties of 
earthquake clustering in time and space. Our main finding is that synthetic catalogs generated by 
the OFC model share practically all properties of real seismicity at a qualitative level, with however 
significant quantitative differences. We find that OFC catalogs can be in large part described by the 
concept of triggered seismicity but the properties of foreshocks depend on the mainshock magnitude, 
in qualitative agreement with the critical earthquake model and in disagreement with simple models 
of triggered seismicity such as the Epidemic Type Aftershock Sequence (ETAS) model g]- Many 
other features of OFC catalogs can be reproduced with the ETAS model with a weaker clustering 
than real seismicity, i.e. for a very small average number of triggered earthquakes of first generation 
per mother-earthquake. Our study also evidences the large biases stemming for the constraints used 
for defining foreshocks and aftershocks. 



I. INTRODUCTION 

Describing and modeling the space-time orga- 
nization of seismicity and understanding the un- 
derlying physical mechanisms to quantify the lim- 
its of predictability remain important open chal- 
lenges. Inspired by statistical regularities such as 
the Gutenberg-Richter Q and the Omori laws, 
a wealth of mechanisms and models have been pro- 
posed. New classes of models inspired or derived 
from statistical physics accompanied and followed 
the proposition, repeated several times under vari- 
ous forms in the last 25 years, that the space-time 
organization of seismicity is similar to the behav- 
ior of systems made of elements interacting at many 
scales that exhibit collective behavior such as in crit- 
ical phase transitions. This led to the concepts of 
the critical earthquake, of self-organized criticality, 
and more generally of the seismogenic crust as a self- 
organized complex system requiring a so-called sys- 
tem approach. 

Our purpose here is to study in depth maybe the 
simplest model of the class of self-organized criti- 
cal models that exhibit a phenomenology resembling 
real seismicity, the so-called OFC sandpile model. 



Real seismicity is usually divided into three ma- 
jor classes of earthquakes, foreshocks, mainshocks, 
and aftershocks. Many different mechanisms have 
been discussed to explain these three classes. The 
OFC model uses only one simple local interaction 
between discrete fault elements but nevertheless ex- 
hibits a sufficiently rich behavior that these three 
classes of events are observed [J]- The motivation 
of our present work is thus to study the main char- 
acteristics of foreshocks and aftershocks in the OFC 
model, and to understand the mechanisms respon- 
sible for earthquake clustering in the OFC model. 
For this, we will interpret our analysis of the OFC 
catalogs in the light of two end-member models 
representing two opposite views of seismicity. The 
first one, mentioned above, is the critical earthquake 
model, which views a mainshock as the special out- 
come of a global self-organized build-up occurring at 
smaller scales. The second end-member is called the 
ETAS (Epidemic- Type Aftershock Sequence) model 
and is nothing but a phenomenological construction 
based on the well-documented empirical Gutenberg- 
Richter law, the Omori law and the aftershock pro- 
ductivity law. In the ETAS model, all earthquakes 
are put on the same footing, any earthquake can 
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trigger other earthquakes and be triggered by other 
earthquakes and the distinction between foreshocks, 
mainshocks and aftershocks is only for convenience 
and does not reflect any genuine difference. The 
ETAS model can thus be considered as a realistic 
statistical null-hypothesis of seismicity. 

The plan of our paper is the following. Section 
m] presents the OFC model, section Hill summarizes 
the phenomenology of real seismicity, and section lTVI 
presents the critical earthquake model and the ETAS 
model. Section M presents the results obtained for 
the OFC model, that are compared with real seis- 
micity and with the two reference models (ETAS and 
critical earthquake models) . We discuss in section lVll 
possible mechanisms for foreshocks and aftershocks 
in the OFC model. Section IVlTI concludes. 

II. THE OLAMI-FEDER-CHRISTENSEN 
(OFC) MODEL 

The Olami-Feder-Christensen (OFC) model is 
defined on a discrete system of blocks or of fault 
elements on a square lattice, each carrying a force. 
The force Fi of a given element i that exceeds a fixed 
threshold Fc (taken equal to 1 without loss of gen- 
erality) relaxes to zero. Such a toppling increments 
the forces on its nearest neighbors by a pulse which 
is a {a < 1/4) times the force Fi of the unstable 
element: 

F^>F^^{f +aF « 

This loading to nearest-neighbors can in turn desta- 
bilize these sites, creating an avalanche. Between 
events, all F^'s increase at the same constant rate, 
mimicking a uniform tectonic loading. The OFC 
model can be obtained as a sandpile analogue of 
block-spring models with the interpretation 



where rii is the actual number of neighbors of site 
i. n.i is always 4 in the case of a square lattice 
with rigid-frame boundaries. For free boundary con- 
ditions used in this study = 4 in the bulk and 
rii = 3 (resp. 2) at the boundaries (resp. corners). 
The symbol k denotes the elastic constant of the up- 
per leaf springs, measured relatively to that of the 
other springs between blocks. The OFC model is 
conservative for fc = for which a = 0.25 and is 
non-conservative for fc > for which a < 0.25. In 
the following, we will compare results obtained for 



fc = 0.5,1,2 and 4, that is for a = 0.222,0.2,0.167 
and 0.125. 

With open or rigid boundary conditions, this 
model seems to show self-organized criticality (SOC) 
^, 8, 9] even in the non-conservative case a < 0.25. 
SOC is the spontaneous convergence of the dynam- 
ics to a statistically stationary state characterized 
by a time-independent power law distribution of 
avalanches. The size of an avalanche is taken to be 
the spanned area s. The underlying mechanism for 
SOC seems to be the invasion of the interior by a 
region spreading from the boundaries, self-organized 
by the synchronization or phase-locking forces be- 
tween the individual elements |lOj |. 

Long-term correlation between large events have 
been documented but, only very recently, Ref. Q 
found the occurrence of genuine sequences of fore- 
shocks and aftershocks that bear similarities with 
real earthquake catalogs. This discovery is interest- 
ing because it suggests that a unique mechanism is 
sufficient to produce a Gutenberg-Richter-like distri- 
bution as well as realistically looking foreshocks and 
aftershocks, without the need for viscous crust re- 
laxation or other mechanisms. Similarly, Ref. [llj | 
found critical precursory activity and aftershock se- 
quences in a sandpile model. However, the precur- 
sors and aftershocks resulted from the interplay be- 
tween the built-in hierarchy of domains and a conser- 
vative sandpile dynamics. The remarkable observa- 
tion of Ref. jjy is that such a hierarchy is not needed 
for foreshocks and aftershocks to occur, when the 
sandpile dynamics is dissipative. However, a hierar- 
chy of faults may be needed to obtain a larger num- 
ber of triggered events than found for the OFC model 
which is more compatible with real seismicity. Our 
goal here is to investigate in details the properties of 
the foreshocks and aftershocks in such a situation, 
that is, in the OFC model. 

Our simulations presented below are performed 
in two-dimensional square lattices L x L with L = 
512, 1024 and 2048. Let us give a correspondence 
between time and space units in the OFC model and 
in the real seismicity. If we consider that the lattice 
of L = 1024 blocks represents a fault of 100 x 100 
km (we neglect the asymmetrical aspect ratio of real 
faults), which is able to produce an earthquake with 
magnitude M — 8, the minimum earthquake of size 
s = 1 generated by the OFC model has a length of 
w 0.1 km, corresponding to a magnitude 2 earth- 
quake. Ref. Q documents that the typical wait- 
ing time between two events involving at least 1000 
blocks is 10~^ in the time unit of the model. This 
event of size s = 1000 corresponds to a rupture 
length of « 3 km or a magnitude Af « 5 earthquake. 
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There are about three events per year with magni- 
tude M > 5 in Southern Cahfornia. This gives the 
correspondence between time units timeoFC in the 
OFC model and in real seismicity: 

timeoFC = 10~^ timCrcai — 100 days . (3) 

Ref. Q describe a particularly active aftershock se- 
quence following an event moving 256^ blocks (mag- 
nitude « 7) lasting timeoFC — 4 • 10^® containing 
about 2500 events. According to Q, this corre- 
sponds to an aftershock sequence lasting only one 
hour. The OFC model clearly does not have the 
time and space scales right to describe the seismic- 
ity in California or in other regions of the world. 
The number of aftershocks is much smaller in the 
OFC model than in California, and aftershocks oc- 
cur at very short times after the mainshock. In the 
OFC model, aftershocks start at topc = 10^ after 
the mainshock (resulting from the numerical preci- 
sion of the simulation). In real seismicity, aftershocks 
are only observed after a few minutes following large 
M = 7 mainshocks, due to the duration of the main- 
shock rupture and to the saturation of the seismic 
network. In real seismicity, aftershocks are also ob- 
served over much longer times, from months to years. 
The aftershock sequences in the OFC model resem- 
ble perhaps more those of deep earthquakes, which 
have few aftershocks ^2 ^| • 



III. PHENOMENOLOGY OF REAL 
SEISMICITY 

The empirical properties of real seismicity dis- 
cussed in this paper are the following. 

1. The Gutenberg- Richter (GR) law [H states 
that the density distribution function P{m) of 
earthquake magnitudes m is 

Po{m)dm = b In(lO) IQ-^f™-™'') dm , (4) 

with a &- value usually close to 1. irid is usually 
a lower bound magnitude of detection, such 
that Pq (m) is normalized to 1 by summing over 
all magnitudes above m^. 

The qualifying property of SOC in the OFC 
model is the existence of a GR-like distribution 
of avalanches sizes. There are several measures 
of sizes. If we take the size defined as the area 
s spanned by an avalanche, the distribution of 
event sizes is exactly given by (0}, where the 
magnitude m is defined by 



The number of topplings is not exactly equal 
to s since a site can topple more than once in a 
given avalanche, possibly being reloaded dur- 
ing its development. However, the difference is 
negligible for our purpose. No multiple relax- 
ations were observed for fc > 2. For k ~ 1, less 
than 1 multiple relaxation per 100,000 earth- 
quakes was found. 



2. The (modified) Omori law [J,llJ,ll5j describes 
the decay of the seismicity rate triggered by a 
mainshock with the time t since the time tc of 
the mainshock 



Na{t) 



K„ 



(6) 



with an Omori exponent pa close to 1. This de- 
cay law can be detected over time scales span- 
ning from weeks up to decades. The time shift 
constant c ensures a finite seismicity rate just 
after the mainshock and is often of the order 
of minutes. 

3. The inverse Omori law 0, ITgI It^ describes 
the average increase of seismicity observed be- 
fore a mainshock and is like © with t — tc 
replaced by ^ 



Nf{t) = 



Kf 



{tc-t + c)Pf ' 



(7) 



with the inverse Omori exponent vt usually 
close to or slightly smaller than pa [13 ■ 

In contrast with the direct Omori law which 
can be clearly observed in a single sequence fol- 
lowing a mainshock, the inverse Omori law can 
only be found by stacking many foreshock se- 
quences because there are huge fluctuations of 
the rate of foreshocks for individual sequences 
preventing the detection of any acceleration of 
seismicity. The well-defined acceleration {TJ of 
the seismicity preceding mainshocks emerges 
only when averaging over many foreshock se- 
quences psj . 

4. The productivity law documents the growth of 
the number of triggered aftershocks as a func- 
tion of the magnitude m of the mainshock: 



Nairn) 



10" 



(8) 



m = logio(s) 



(5) 



where the exponent aa is usually found in the 
range 0.7—1 (see and references therein). 
This value of the exponent aa may reveal a 
fractal spatial distribution of aftershocks fisll. 
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5. Aftershock diffusion. Several studies have re- 
ported "aftershock diffusion," the phenomenon 
of expansion or migration of aftershock zone 
with time [20ll2ll |. However, the present state 
of knowledge on aftershock diffusion is confus- 
ing because contradictory results have been ob- 
tained, some showing almost systematic diffu- 
sion whatever the tectonic setting and in many 
areas in the world, while others do not find ev- 
idences for aftershock diffusion [23, ■ The 
shift in time from the dominance of the after- 
shock activity clustered around the mainshock 
at short times after the mainshock to the delo- 
calized background activity at large times may 
give rise to an apparent diffusion of the seis- 
micity rate when using standard quantifiers of 
diffusion processes [23 ■ 

6. Foreshock migration. Foreshock migration to- 
wards the mainshock as time increases up to 
the time of the mamshock has also been doc- 
umented (3 13 101 but may be due to an 
artifact of the background activity, which dom- 
inates the catalog at long times and distances 
from the mainshock p5l |. Indeed, by an argu- 
ment symmetrical to that for aftershocks, the 
shift in time from the dominance of the back- 
ground activity at large times before the main- 
shock to that of the foreshock activity clustered 
around the mainshock at times just before it 
may be taken as an apparent inverse diffusion 
of the seismicity rate. 

7. The average distance Ra between aftershocks 
and the mainshock rupture epicenter has been 
found to be proportional to the rupture size of 
the mainshock, leading to the following scaling 
law 



Ra - 10' 



a. 5m 



Q.5 



(9) 



relating Ra and the mainshock magnitude m or 
the mainshock rupture surface s. A similar law 
is believed to hold for the average distance _Rf 
between foreshocks and the mainshock [271 l28l| 

8. Foreshock magnitude distribution. Many stud- 
ies have found that the apparent 6- value of the 
magnitude distribution of foreshocks is smaller 
than that of the magnitude distribution of the 
background seismicity and of aftershocks (see 
[l^ and references therein). 

9. Number of foreshocks and aftershocks per 
mainshock. Foreshocks are less frequent than 
aftershocks Q 0, [^^ • The ratio of foreshock 



to aftershock numbers is in the range 2-4 for 
M = 5 — 7 mainshocks, when selecting fore- 
shocks and aftershocks at a distance in the 
range R = 50—500 km from the mainshock and 
for a time in the range T — 10— 100 day s before 
or after the mainshock [llllllllillli. 



IV. END-MEMBER MODELS OF 
SEISMICITY: ETAS AND CRITICAL 
EARTHQUAKE MODELS 

A. The ETAS model 

The ETAS (Epidemic- Type Aftershock Sequence) 
model was introduced in 2] and in (SOj (in a slightly 
different form). Contrary to what its name may im- 
ply, it is not only a model of aftershocks but a general 
model of seismicity. 

This model of multiple cascades of earthquake 
triggering avoids the division between foreshocks, 
mainshocks and aftershocks because it uses the same 
laws to describe all earthquakes. Because of its sim- 
plicity, it is natural to consider it as a null hypothesis 
to explain the OFC catalogs and real data. Its choice 
as a reference is also natural because it is a noth- 
ing but a branching model of earthquake interac- 
tions and can thus be considered as a mean field ap- 
proximation of more complex interaction processes. 
Branching processes can also be considered as natu- 
ral mean field approximations of SOC models and in 
particular of the OFC model jsj (see also chapter 
15 in 0). The approximation consists usually in the 
fact that branching models neglect the loops occur- 
ring in the cascade characterizing a given avalanche. 
Note that standard branching models study the de- 
velopment of a single avalanche while ETAS describe 
a catalog of earthquakes. 

The ETAS model uses three of the above empiri- 
cal laws as direct inputs (Gutenberg-Richter law 
Omori's law © and aftershock productivity law (jSJ). 
In the ETAS model, a main event of magnitude m 
triggers its own primary aftershocks (considered as 
point processes) according to the following distribu- 
tion in time and space 



0™(r,i) dr dt = K 10° 



e dt 



df" dr 



{t + c)^+o (r + rf)i+^ ' 
(10) 

where r is the spatial distance to the main event. 
The spatial regularization distance d accounts for the 
finite rupture size. The power law kernel in space 
with exponent ^ quantifies the fact that the distri- 
bution of distances between pairs of events is well 
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described by a power-law ^3^. The ETAS model 
assumes that each primary aftershock may trigger 
its own aftershocks (secondary events) according to 
the same law, the secondary aftershocks themselves 
may trigger tertiary aftershocks and so on, creat- 
ing a cascade process. The exponent 1 + 6* is not 
the observable Omori exponent pa but defines the 
"bare" Omori law for the aftershocks of first genera- 
tion. The whole series of aftershocks, integrated over 
the whole space, can be shown to lead to a "renor- 
malized" (or dressed) Omori law, which is the total 
observable Omori law '33^. To prevent the process 
from dying out, a small Poisson rate of uncorrelated 
seismicity driven by plate tectonics is added to rep- 
resent the effect of the tectonic loading in earthquake 
nucleation. 

The ETAS model predicts the following properties. 

1. The same productivity law © is found for the 
total number of events (and not just from the 
first generation) ,33]. This law describes an 
average power law increase of the number of 
directly and indirectly triggered earthquakes 
per triggering earthquake as a function of the 
mainshock size s. 

2. The ETAS model finds a global direct Omori 
law (for aftershocks) different from the "bare" 
Omori law - l/{t + c)i+^ defined in ^ for 
the first generation aftershocks, with a renor- 
malized exponent Pa smaller than from 1 + 9 

3. An inverse Omori law for foreshocks is found to 
result simply from the existence of the direct 
"bare" Omori law 110|l and from cascades of 
multiple triggering 18J . 

4. The ETAS model predicts a modification of the 
magnitude distribution before a mainshock 
of magnitude mj\/ , characterized by an increase 
of the proportion of large earthquakes accord- 
ing to the following expression [13 : 

P{m\mM) = (1 - Q{t))Po{m) + Q{t)dP{m) , (11) 

where Pq (™) is the unconditional GR distribu- 
tion 

dP{m) - lQ-b'ira-r,ia) ^ ^i^J^ b' = b - Ua , (12) 

and 



tc is the time of the mainshock and C is a 
numerical constant. The prediction (|ll|l with 
(|12I13|I is that the magnitude distribution is 
modified upon the approach of a mainshock 
by developing a bump in its tail which takes 
the form of a growing additive power law con- 
tribution with a new 6'-value. This prediction 
has been validated very clearly by numerical 
simulations of the ETAS model ^3 ■ 

5. In the ETAS model, the properties of fore- 
shocks are independent of the mainshock mag- 
nitude, because the magnitude of each event 
is not predictable but is given by the GR law 
with a constant 6- value |lq| . 

6. By the mechanism of cascades of triggering, 
the ETAS model also predicts the possibility 
for large distance and long-time build-up of 
foreshock activity as well as the migration of 
foreshocks toward mainshocks. [2^ l25l | . 

These predictions of the ETAS model are in good 
agreement with observations of the seismicity in 
Southern California psf . 

B. The critical earthquake model (CEM) 

Maybe the first work on accelerated seismicity 
leading to the concept of criticality is .27], who ob- 
served that the trailing total sum of the source areas 
of medium size earthquakes accelerates with time on 
the approach to a large or great earthquake. The 
theoretical ancestor of the critical earthquake con- 
cept can probably be traced back to 0], who used 
a branching model to illustrate a cascade of earth- 
quake ruptures culminating in complete collapse in- 
terpreted as a great one. Ref. (35i] proposed a renor- 
malization group analysis of a percolation model of 
damage/rupture prior to an earthquake paralleling 
p6>] . which emphasized the critical point nature of 
earthquake rupture following an inverse cascade from 
small to large scales. Refs. [37ll38| | are probably the 
first ones to introduce the idea of a time-to-failure 
analysis in the form of a second order nonlinear dif- 
ferential equation, which for certain values of the pa- 
rameters leads to a solution of the form of a time- 
to-failure equation describing the power law acceler- 
ation of an observable with time: 

eit)ocA-B{t,~ty , (14) 



where e{t) is for instance the cumulative Benioff 
strain (square root of earthquake energy), A and 
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B are positive constants, is the critical time of 
the mainshock and < z < 1 is a critical ex- 
ponent. Ref. 39] introduced equation (|14|l to fit 
and predict large earthquakes. Their justification 
of p4|l was a mechanical model of material dam- 
age. Ref. did not mention the critical earth- 
quake concept. Ref. ji^] proposed to reinterpret 
[39| and all these previous works and to gener- 
alize them using a statistical physics framework. 
The concept of a critical earthquake described in 
Ref. ^3 corresponds to viewing a major or great 
earthquake as a genuine critical point in the statis- 
tical physics sense. In a nutshell, a critical point 
is characterized by long-range correlations and by 
power laws describing the behavior of various ob- 
servables on the approach to the critical point. This 
concept has been elaborated in subsequent works 
[HISIillilEllillillillillii The criti- 
cal earthquake model predicts a power-law increase 
of the number, of the cumulative displacement and of 
the average energy of foreshocks before large earth- 
quakes. According to this model, the modification 
of seismic activity should be more apparent before 
larger mainshocks. Therefore, we should measure a 
positive value for a/ characterizing foreshocks. The 
critical model also predicts that the preparation zone 
or foreshocks cluster size i?/ should increase with 
the mainshock size s as Rf s'^f with qf > 0, as 
observed by [27i (see ^3 f^'' extended compi- 
lation and discussion). 

V. PROPERTIES OF THE SYNTHETIC 
SEISMICITY GENERATED BY THE OFC 
MODEL 

A. Definitions of foreshocks and aftershocks 

The first striking observation, extending the dis- 
covery of , is that all the properties of real seismic- 
ity discussed in section Hill are found to exist at least 
at a qualitative level in synthetic catalogs generated 
by the OFC model, as we now show. 

Deciding what is a foreshock or an aftershock is 
not straightforward, and several methods have been 
proposed to define aftershocks or foreshocks [23L lilj . 
A clear identification of foreshocks, aftershocks and 
mainshocks is hindered by the fact that nothing dis- 
tinguishes them in their seismic signatures: at the 
present level of resolution of seismic inversions, they 
are found to have the same double-couple structure 
and the same radiation patterns ,51j . It is thus im- 
portant to consider several alternative definitions. 
• Definition d — 0. The common definition of 



foreshocks and aftershocks is as follows. A main- 
shock is first defined as the largest event in a se- 
quence. Foreshocks and aftershocks are then defined 
as earthquakes in a pre-specified space-time domain 
around the mainshock, and are thus constrained to 
be smaller than the mainshock. The space-time win- 
dow used for foreshock and aftershock selection is 
defined such that the seismicity rate is much larger 
than the background seismicity. For simplicity, we 
choose a space window that covers the whole lattice, 
and we choose a time window with a duration equal 
to 1% of the recurrence time of the mainshock, so 
that the seismicity within this window is much larger 
than the background seismicity. 

However, recent empirical and theoretical studies 
suggest that this definition might be arbitrary and 
physically artificial [H [H HI US [H HI In- 
deed, the magnitude of an earthquake seems to be 
unpredictable |l9j . therefore the same mechanisms 
responsible for the triggering of small earthquakes 
(usual "aftershocks" for d = 0) may also explain the 
triggering of larger earthquakes (defined as "main- 
shocks" for d = 0). This hypothesis has been tested 
in (25j and accounts very well for the properties of 
earthquake clustering in real seismicity. We thus use 
two other definitions of foreshocks and aftershocks, 
which do not constrain aftershocks and foreshocks to 
be smaller than the mainshock, in order to test how 
the selection procedure impacts on foreshock and af- 
tershock properties. 

• Definition d — 1. Mainshocks are now defined as 
any earthquake that was not preceded by a larger 
earthquake within a time window equal to 1% of the 
mainshock recurrence time, but can be followed by a 
larger event. This rule aims to select as aftershocks 
the events that have been triggered directly or in- 
directly by the mainshock, removing the influence 
of large earthquakes that occurred before the main- 
shock. We use the same space-time domain as for 
d = to select foreshocks and aftershocks around 
the mainshock. Foreshocks are thus smaller than 
mainshocks but aftershocks can be larger. 

• Definition d = 2. Same as d = 1 without any 
constraint on the aftershocks and mainshocks: each 
event of appropriate size is a mainshock, preceded by 
foreshocks and aftershocks which are selected solely 
on the basis of their belonging to the pre-specified 
space-time domain around the mainshock, with ar- 
bitrary magnitudes. For foreshocks, this corresponds 
to the "Type 11" foreshocks introduced in 0, in 
order to remove any spurious dependence of fore- 
shock properties on the mainshock magnitude. 

For d = 2, each earthquake can be a foreshock or 
an aftershock to several mainshocks and can also be 
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a mainshock. 



B. Direct (aftershocks) and inverse 
(foreshocks) Omori laws 

Figure ^ shows 15 individual sequences of fore- 
shocks (bottom) and aftershocks (top), for main- 
shocks of size s > 2048 with more than 1000 fore- 
shocks or aftershocks, generated in a system of size 
L — 2048, dissipation index k — 2 and selected with 
definition d = 0. The direct Omori law © is strik- 
ingly clear for each individual sequence while, for 
foreshocks, there is almost no increase of the seis- 
mic rate for individual sequences: the inverse Omori 
law (UI) is only observed when stacking many se- 
quences, like in the ETAS model 0|. This is the 
clearest indication that the foreshock activity may 
be better described by a cascade ETAS-type model 
than by a critical earthquake model (but see below 
for other observations that may modulate this pre- 
liminary conclusion). The OFC model shares this 
property with empirical seismicity and with the 
ETAS model 

We now describe the results obtained by aver- 
aging over a large number of sequences, which al- 
low us to decrease the noise level and to look at 
smaller mainshocks. We have generated synthetic 
catalogs with the OFC model, using a lattice size 
L = 128, 256, 512, 1024, 2048 and different values of 
the dissipation index k = 0.5, 1,2, 4. For each cata- 
log, we have selected aftershocks and foreshocks fol- 
lowing the different definitions d = 0, 1,2 explained 
in section IV Al We have then stacked all sequences 
by superposing all sequences translated in time so 
that the mainshock occurs at time t ^ 0. We have 
first analyzed the change of the seismicity rate be- 
fore and after a mainshock. For each range of main- 
shock size s between 2 and 2^^, increasing by a fac- 
tor 2 between each class, we compute the average 
seismicity rate Nf{t) and Na{t) as a function of the 
time before and after the mainshock. The results for 
L = 2048, k = 2 and d = are illustrated in Figure 
121 The rate of aftershocks obeys the direct Omori 
law (0 and the increase of the seismicity rate ob- 
served when averaging over many sequences follows 
the inverse Omori law {TJ. We measure the Omori 
exponents pf and pa by fitting the number Na{t) 
of aftershocks and the number Nf{t) of foreshocks 
by a power law ~ l/|t|P"'^ using a linear regression 
of In(A^) as a function of ln|t|, in the time interval 
\t\ > 5- 10~^^ and \t\ < tmax, where the upper bound 
tmax is given by the condition that the seismicity rate 
is much larger than the background level. 



The duration of aftershock or foreshock sequences 
is much smaller than the average recurrence time 
of earthquakes. This is also shown in Figure 1 of 
Ref. P|, which represents the distribution of inter- 
event times. A deviation from the Poisson distribu- 
tion is evident only for small times compared to the 
average inter-event times. 

Table provides the values of the exponents Pa 
and Pf as a, function of L, k and d. We find similar 
Omori exponents for foreshocks and aftershocks with 
Pf < Pa < 1- We observe the same time-dependence 
of the seismicity rate (same exponents Pa and pf) 
for all mainshock sizes, only the absolute value of 
the seismicity rate depends on the mainshock magni- 
tude. The exponents pa and pf are found to increase 
with k from pa ~ 0.5 for k = 0.5 to Pa ~ 0.9 for 
k — 4, but the duration of the aftershock and fore- 
shock sequence does not change significantly with k. 
The Omori exponents do not depend on the rules 
of selection d. The number of foreshocks and after- 
shocks thus increases if the dissipation increases, and 
is almost negligible in the non-dissipative case. 

C. Dependence of the number of aftershocks 
and foreshocks with the mainshock size 

Figure 13 represents the dependence of the number 
of aftershocks and foreshocks with the mainshock 
size, for different rules of selection d = 0,1,2. We 
observe a power-law increase of the number of fore- 
shocks Kf and aftershocks Ka defined in 10 and 
with the mainshock size s according to 

Ka{s) - s"" , and Kf{s) ~ s°'f , (15) 

for s < 10'^. The exponents aa and af measured for 
s < 1000 increase with the dissipation index k (see 
Table 1). The results are very similar for d — 1 and 
d = 0. The exponent aa is slightly smaller for d — 1 
than for d = 0, because for d = we impose the 
aftershocks to be smaller than mainshocks. Small 
events are more likely to trigger an event larger than 
themselves than larger mainshocks and thus to be 
rejected from the analysis. Therefore the rule d = 
underestimates the number of earthquakes triggered 
by small mainshocks. 

For d — 2 and for small s, Ka{s) and Kf{s) are 
much larger than with d = 0, and are almost inde- 
pendent of s for s < 100. This results from the fact 
that, for d — 2 a significant fraction of "mainshocks" 
are triggered by a previous larger event, and thus 
the events classified as aftershocks may be in fact 
triggered by earthquakes that occurred before their 
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"mainshock" . The results obtained with d = 2 re- 
covers those obtained with d = and d = 1 for large 
s. The correct value of the exponent aa of the after- 
shocks productivity law is thus the value obtained 
for d = 1. 

The number of foreshocks is generally smaller than 
the number of aftershocks and increases more slowly 
with s (a/ < Ua)- 

For large mainshock sizes s > 1000, we observe a 
saturation of Ka and Kf, and the numbers of fore- 
shocks and aftershocks increase slower with s than 
predicted by (|15|l . This saturation size does not de- 
pend neither on k nor on L. The effect of the system 
size L is only to change the shape of the functional 
form of Ka{s) and Kf{s) for large s: the saturation 
of Ka.f{s) for s > 1000 is more obvious for smaller 
L. 



D. Spatial distribution of foreshocks and 
aftershocks 

Figure ^ shows the stress field immediately be- 
fore and after a large mainshock. Following the 
mainshock, many elements on the boundary of the 
avalanche and within the avalanche have been loaded 
by the rupture, and are likely to generate aftershocks 
after the mainshock. There are a few large patches 
of elements within the avalanche that did not break 
during the mainshock, as illustrated in the lower 
panels of Figure 0] but most aftershocks initiate on 
smaller clusters of a few unbroken elements shown 
as white spots in the central lower panel of Figure 
01 The density of such white spots observed in this 
square is typical of the rest of the stress field over 
the area spanned by the mainshock. In the language 
of seismicity, such white spots are "asperities" which 
carry a large stress after the mainshock and are nu- 
cleation point for future aftershocks. These white 
spots are also found on the avalanche boundaries. 

The aftershock cluster size Ra, defined as the av- 
erage distance between the mainshock and its af- 
tershocks, is close to the mainshock size i = ^/s 
at small times, and then increases to a value close 
to the system size L at large times (see Figure (SJ. 
This crossover of Ra{t) corresponds to a transition 
between the aftershock activity at small times and 
short distances from (or within the area spanned 
by) the mainshock, to the uncorrelated seismicity at 
larger times after the mainshock and with a uniform 
spatial distribution. We observe the same pattern 
for the spatial distribution of foreshocks. The time 
of the crossover of Ra{t) and Rf{t) increases with 
s because the number of foreshocks and aftershocks 



increases with s and thus remain for a longer time 
above the uncorrelated background seismicity. In the 
short times regime where uncorrelated seismicity is 
negligible, we find a very weak, if any, diffusion of 
aftershocks, as measured by an increase Ra ~ t^" of 
the aftershock zone size with the time after the main- 
shock. Similarly, we observe a very weak migration 
of foreshocks toward the mainshock, characterized by 
a decrease i?/ ~ (to — t)^f of the foreshock zone as 
the mainshock approaches. The Hurst exponents Ha 
and Hf are not statistically different from 0, showing 
that the size of the foreshock and aftershock zones 
do not change significantly with time. 

We observe on Figure a power-law increase of 
aftershock zone size and of the foreshock zone size 
with the mainshock size according to 

-Ra(s)~s«", and i?/(s)-s«^, (16) 

in the range 10 < s < 10*. The exponents qa and 
Qf are given in Table 1 as a function of k, L, and d. 
In contrast with real seismicity, the aftershock zone 
area is not proportional to the mainshock size s"'^, 
but it increases slower with s {qa ~ 0.3 < 0.5 for 
d = or 1). This is probably due to the effect of sec- 
ondary aftershocks, which increase the effective size 
of the aftershock zone for small mainshocks. Sec- 
ondary aftershocks are more important for d = 1 
than for d — 0, which explains why qa is smaller for 
d = 1 than for d = 0. The average value of the 
foreshock zone ( "zone of mainshock preparation" ) is 
smaller than the aftershock zone, except for defini- 
tion d — 2. 

An increase of i?/ with s according to H16|) has 
been reported by Ref. [23, for individual se- 
quences, with an exponent (?/ = 0.44 28], but was 
not observed by Ref. '25'| when using stacked se- 
quences and when allowing foreshocks to be larger 
that the mainshock {d — 2). This increase of i?/ 
with s is not observed by the ETAS model, because 
the magnitude of each earthquake is drawn at ran- 
dom, independently of previous seismicity, and thus 
all properties of foreshocks must be independent of 
the mainshock size. 

For d = 2, the average zone sizes Ra and of i?/ are 
much larger than for d = 0, because we include fore- 
shocks and aftershocks larger than the mainshock. 
^For s < 1000, the values of Ra and Rf are almost 
constant of the order of Raj ~ 100. This may reflect 
the fact that, for d — 2 and for small mainshocks, 
most aftershocks are not triggered by the mainshock 
but by a previous larger event. 
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E. Distribution of avalanche sizes 

For the whole catalogs, the distribution of event 
sizes is a power-law characterized by an exponent 
b, with an exponential roll-off for large sizes s > 
10~^L^. Table^shows that the exponent b increases 
as the dissipation index increases, from b « 0.7 for 
k = 0.5 to the realistic value b « 0.95 for k = A, in 
agreement with previous works 0, |^ . 

Figures [7| test the stationarity of the magnitude 
distribution for foreshocks and aftershocks. The 
magnitude distributions for foreshocks and after- 
shocks have been fitted with expression Hll|l and 
these fits are shown as the colored lines. The devia- 
tion of the magnitude distribution from the average 
GR law for foreshocks and aftershocks are well rep- 
resented by Ullfl with b' in the range 0.3 — 0.5 and 
with Q{t) increasing as a power law according to ifT^ 
with 1/ « 0.1 for foreshocks and for aftershocks. 

For aftershocks, the time-dependence of Q{t) de- 
scribes a decrease of the deviation from the GR law 
for latter aftershocks. While these fits are good, 
there is an important caveat: the prediction b' = 
6 — a of the ETAS model is not valid here, as the 
OFC model gives aa ~ 0.8 for k = 2, which is 
slightly larger than b — 0.78. This regime aa > b 
gives an explosive behavior in the ETAS model with 
a stochastic finite time singularity 56] unless the GR 
law is truncated by an upper bound or is rolling over 
to a faster decay for the large earthquake sizes. 

The modification of the magnitude distribution for 
aftershocks in the OFC model, is much weaker than 
for foreshocks, but is significant. This implies that 
the magnitude distribution in the OFC model is not 
stationary because the magnitude of triggered earth- 
quakes is correlated with the mainshock magnitude, 
in contradiction with a crucial hypothesis of the 
ETAS model and with real catalogs. Figure |S1 shows 
that, for d = 1 (no constraint on aftershock mag- 
nitudes), the change of the magnitude distribution 
is almost independent on the mainshock magnitude 
(like in the ETAS model). However, there is a larger 
proportion of medium-size events for smaller main- 
shocks than for larger ones. This means that smaller 
mainshocks have a tendency to trigger smaller after- 
shocks than larger mainshocks. 

This result is in contradiction with the ETAS 
model, which does not reproduce a dependence of 
the aftershock magnitude as a function of the size of 
the triggering earthquake or as a function of the time 
since the mainshock. Observations of real seismicity 
[T^ do not show any dependence between the after- 
shock magnitude and the mainshock magnitude, but 



the catalogs available are much smaller than for our 
OFC simulations. 



VI. MECHANISMS FOR FORESHOCKS 
AND AFTERSHOCKS IN THE OFC MODEL 

Is the increase of the number of aftershocks and 
foreshocks with the mainshock magnitude real or is 
it just the result of a selection bias introduced by the 
standard definition d = 0, which requires that main- 
shocks are the largest events in the cluster? Indeed, 
in the OFC catalogs, the clear power-law increase of 
the number of aftershocks and foreshocks with the 
size s of the mainshock, found when we define the 
mainshock as the largest event (d = 0), almost dis- 
appears when aftershocks or foreshocks are not con- 
strained to be smaller than the mainshock {d — 2), 
as shown in Figure 13 The question of the impact of 
the definition is thus essential. 

For foreshocks, we consider two possible interpre- 
tations, the ETAS model described in section HV Al 
and the critical earthquake model (CEM) summa- 
rized in section llV Bl The fact that, using d = 0, 
both the number of foreshocks and the average fore- 
shock cluster size increase with the mainshock size 
s seems to favor the critical model, but Ref. 
have shown that the constraint that foreshocks must 
be smaller than the mainshock (d = 0) leads to an 
artificial increase of the number of foreshocks and 
of the foreshock cluster size with s (see Fig. 5 of 
|0|). However, this spurious increase of the number 
of foreshocks with the mainshock magnitude should 
be observed only for small s, and should not exist 
in the case d = 2 (without constraints on foreshock 
and mainshock magnitudes). The fact that a weak 
increase of the number of foreshocks with s is ob- 
served even for very large s and for all definitions 
d = 0, 1, 2 suggests that the effect is genuine. Such a 
dependence of foreshock properties on the mainshock 
size cannot be reproduced with the ETAS model, but 
suggest that the critical model provides a more rele- 
vant description of these observations. 

The case d = 2 destroys the dependence of the af- 
tershock number with s for small s, which is a real 
physical property of aftershocks because it is also 
observed for d = 1. The same effect may also be at 
work for foreshocks and explain why, if the proper- 
ties of foreshocks are physically dependent on the 
mainshock magnitude, this dependence is not ob- 
served when using d = 2. Definition d = 2 may 
be mixing "critical foreshocks", i.e. events that be- 
long to the preparation phase of the future main- 
shock, with "triggering-triggered" pairs (which are 
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the usual "foreshock-mainshock" pairs in the ETAS 
model). Thus, the absence of conditioning {d = 2) 
seems to destroy the dependence of both the fore- 
shock and aftershock properties with the mainshock 
magnitudes, i.e., aaj ~ and qaj ~ both for 
foreshocks and aftershocks. But, for aftershocks, we 
know for a truth that the dependence of the number 
(ua) and the cluster size (qa) are real and physical 
because they are observed in the case d = I which 
do not constrain aftershocks to be smaller than the 
mainshocks. 

For aftershocks selected according to d = 2 and 
for small mainshock sizes s, the scenario, accord- 
ing to which most aftershocks are not triggered by 
the mainshock per se but by a previous larger event, 
seems to explain both the fact that (i) the number of 
aftershocks is almost constant with s for s < 100 and 
(ii) the size of the aftershock cluster is almost con- 
stant with s for small s. It is however surprising, if 
this interpretation is correct, that we observe a pure 
Omori law for aftershocks and foreshocks in the case 
d — 2, without any change in the Omori exponent 
with s and without any roll-off at small times. 

The observations that earthquake clustering in 
the OFC model is much weaker than for real seis- 
micity implies that the "branching ratio" (average 
number of triggered earthquakes of first generation 
per mot her- earthquake) must be close to for the 
ETAS model to match the OFC model. But such a 
small branching ratio would lead to Omori's expo- 
nent larger than 1, in contradiction with our simula- 
tions. Thus, the pure branching model of triggered 
seismicity is insufficient to fully account for the quan- 
titative aspect of the OFC seismicity. 

In order to better understand the mechanism re- 
sponsible for aftershocks in the OFC model, we have 
imposed several types of perturbations to the nor- 
mal course of an OFC simulation to obtain the re- 
sponse of the system. First, we have simulated ran- 
dom isolated disturbances consisting in choosing ran- 
domly and independently 1024 sites and adding to 
them random amounts of stress drawn in the inter- 
val [0,0.01] or [0,0.1]. Repeating such disturbance 
100000 times, we find no observable seismicity trig- 
gered by this perturbation. This shows that the 
aftershocks require a coherent spatial organization 
over a broad area. In contrast, pasting (or grafting) 
the stress map as shown in figure^lof those sites that 
participated in a given large mainshock in a given 
simulation and their nearest neighbors onto another 
independent simulation gives a perfect Omori law 
following this graft in this second simulation as if the 
foreign stress map of the first simulation was part 
of the second simulation. The resulting aftershock 



sequences are similar to natural sequences. This 
demonstrates that the presence of asperities close 
to the failure threshold inside the boundary of the 
avalanche can produce realistic aftershock sequences 
with an Omori law temporal decay. To investigate 
further if it is the spatial connectivity of the per- 
turbation which is important to get aftershocks, we 
have also raised simultaneously the stress by various 
amounts within squares of size 64 x 64. Performing 
this perturbation 40000 times at random instants, 
we observe that this sometimes results in quite large 
earthquakes immediately following the perturbation, 
but there is no significant activity afterwards, and 
nothing that looks like an Omori law of aftershocks. 
Thus, it seems impossible to generate aftershocks by 
introducing a random or deterministic spatially ex- 
tended perturbation. This suggests that aftershocks 
require not only the occurrence of mainshocks to re- 
distribute stress but also a spatial organization of 
the stress field prior to the mainshock so that many 
"asperities" (as illustrated in figure^ can be created 
by the interplay of the stress organization before the 
mainshock and the stress redistribution by the main- 
shock. 



VII. CONCLUSION 

The results obtained in this paper can be viewed 
from two different perspectives. On the one hand, 
we are adding to the phenomenology of one of the 
most studied model of self-organized criticality. Ex- 
tending Hergarten and Neugebauer's announcement 
0, we have shown evidence that the self-organized 
critical state of the OFC model is much richer than 
previously thought, with important correlation pat- 
terns in space and time between avalanches. We 
have obtained quantitative scaling laws describing 
the spatio-temporal clustering of events in the OFC 
model. On the other hand, we have shown that 
what is probably the simplest possible mechanism 
for the generation of earthquakes (slow tectonic load- 
ing and sudden stress relaxation with local stress 
redistribution) is sufficient to recover essentially all 
known properties documented in seismic catalogs, at 
a qualitative level. Specifically, we have found that 
foreshocks and aftershocks follow the inverse and di- 
rect Omori power laws, as in real seismicity but with 
smaller exponents. We have found the productivity 
power law of aftershocks as in real seismicity. In 
contrast with real seismicity, we also found a power- 
law increase of foreshock productivity with the main- 
shock size. The nucleation of aftershocks at "asper- 
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ities" located on the mainshock rupture plane or on 
the boundary of the avalanche is also in agreement 
with seismic observations. These findings are inter- 
esting because they add on the list of possible mech- 
anisms for foreshocks and aftershocks that have been 
discussed previously in the literature. 

We also found that the predictability in- 
creases with the mainshock magnitude in the 
OFC model, because the number of foreshocks 
seems to increase with the magnitude of the 
mainshock, a feature that is not observed in 
real seismicity 25j. See the debate in Nature at 
I http: / / helix. nature. com/debates/earthquake/ and 
|57l 15^ for related discussions on the predictability 
of real earthquakes and the use of models from 
statistical physics. 

While most of the OFC dynamics can be qual- 
itatively captured by the simple ETAS model of 
triggered seismicity, this property that the number 
of foreshocks seems to increase with the magnitude 
of the mainshock is better explained by the criti- 
cal earthquake model. We have systematically com- 
pared the statistical properties of the avalanches gen- 
erated by the dynamics of the OFC model with those 
predicted by the ETAS model on the one hand and 
by the critical earthquake model on the other hand. 
These two models constitute end-member models of 
seismicity. We can conclude that, while the de- 



tailed quantitative predictions of the ETAS model 
for its application to the OFC model are not exactly 
right, the physical mechanism of cascades of trig- 
gered events inherent in the ETAS model captures 
at least qualitatively the observed regularities found 
in the OFC model. This suggests a picture in which 
future avalanches are triggered by past avalanches 
through "asperities" located either within the plane 
of past avalanches or at their boundaries. In addi- 
tion, we find that this triggering mechanism presents 
a degree of cooperativity, as the number of foreshocks 
increases with the mainshock size. In other words, 
this suggests that asperities interact via avalanches, 
and when their number and size increase in a given 
location, they can produce larger avalanches. 
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FIG. 1: 15 individual sequences of aftershocks (a) and 
foreshocks (b) , for mainshocks of size s > 2048 with more 
than 1000 foreshocks or aftershocks, generated in a sys- 
tem of size L — 2048, dissipation index k — 2 and se- 
lected with definition d = 0. 



TABLE I: Aftershock and foreshock properties in the 
OFC model as a function of the system size L and of the 
dissipation index k, for different definitions of foreshocks 
and aftershocks d (see section IV Al . b is the exponent of 
the cumulative distribution of avalanche sizes, measured 
for the whole catalog and discussed in section IV El pa 
and p/ are the exponents of the direct and inverse Omori 
laws (see section FyBt . Oa and a/ characterize the de- 
pendence of the aftershock and foreshock rate with the 
mainshock size s (see section IV CI . qa and g/ describe 
the scaling of the aftershocks and foreshock zone size 
with s (see section IV Dt . A^max is the maximum rate of 
aftershocks at t = 10^^^ normalized by the background 
rate. 
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FIG. 2: Direct ((a), aftershocks) and inverse ((b), fore- 
shocks) Omori law for synthetic catalogs generated with 
the OFC model with L — 2048, dissipation index k = 2 
and definition d = 0. The seismicity rate is normalized 
by the background rate and by the number of mainshocks 
in each class. The mainshock size increases from s — 2 
(bottom curve) to s — 2^^ (top curve) with a factor 2 
between each curve. 
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FIG. 3; Productivity of aftershocks Ka (circles) and 
foreshocks Kf (crosses) as a function of the mainshock 
size s for synthetic catalogs generated with the OFC 
model with L = 2048, dissipation index k = 2 and defi- 
nition d = (a), d = 1 (b) and d = 2 (c). 
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FIG. 4: Stress field immediately before (a) and after (b) 
a mainshock. The stress change due to the mainshock 
in shown in (c). The elements that broke during the 
avalanche are shown in dark in (c) (stress decrease) , and 
were mostly close to the rupture threshold before the 
mainshock (light gray in (a)). The upper panels show the 
whole grid of size L = 1024 and the lower plots represent 
a subset of the grid delineated by the square in the upper 
plot. 
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FIG. 5: Average distance between mainshocks and after- 
shocks (a) and mainshocks and foreshocks (b) as a func- 
tion of the time after (a) or before (b) the mainshock, for 
synthetic catalogs generated with the OFC model with 
L = 2048, dissipation index k — 2 and definition d = 0. 
Each curve corresponds to a different mainshock size size 
increasing from s — 2 (bottom curve) to s = 2^^ (top 
curve) with a factor 2 between each curve. 
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FIG. 6: Average distances Ra (circles) and Rj (crosses) 
as a function of the mainshock size s, for k — 2, L — 2048, 
and for definitions d = (a), d = 1 (b) and d = 2 (c). 
The continuous lines are fits of Ra,f{s) by a power-law 
s". 




FIG. 7: Aftershock (a) and foreshock (c) size distribu- 
tion as a function of the time before or after the main- 
shock, for k — 2 and L — 2048, and for a mainshock size 
in the range s = 2048 — 4096. Aftershocks were selected 
with d — 1 and foreshocks with d = 2 (without constrains 
on the magnitude of foreshocks and aftershocks). The 
avalanche size distributions are constructed with a loga- 
rithmic binning (linear bin in magnitudes), whose slope 
gives the GR 6-value. The black line in (a) in (c) shows 
the size distribution for the whole catalog for reference. 
The different colors from blue to red correspond to dif- 
ferent time windows closer and closer to the mainshock, 
from \tc-t\ = 5x10"^ to \tc-t\ = 5x10"^". The colored 
lines in (a) and (c) are the fit of the foreshock and after- 
shock magnitude distribution P{jn) (with m — logjQ(s) 
as defined in @) by expression Ijllfl . The corresponding 
values of b'{t) il2l and Q{t) 11311 are shown in (b) and 
(d) for foreshocks (crosses) and aftershocks (circles), as 
a function of time before/after the mainshock. 
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FIG. 8: Aftershock (a) and foreshock (b) size distri- 
bution for A: = 2 and L = 2048. Aftershocks [d = 1) 
and foreshocks (d = 2) are selected in a time window 
5 X lO"^'^ < |ic — i| < 5 X 10~^°. Each curve corresponds 
to a different mainshock size between s = 4 (blue curve) 
and s — 2^^ (red curve) increasing by a factor 4 between 
each curve. The black line shows the (unconditional) 
distribution of the whole catalog. 



